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Abstract. Fluctuations of cell state, e.g., abundances of some proteins, have 
attracted much attention both theoretically and experimentally. The distribution of 
such state over cells, however, is not only a result of intracellular stochastic process, 
but is also influenced by the growth in cell numbers that depends on the state. By 
incorporating the growth-death process into the standard Fokker-Planck equation for 
the probability distribution, a nonlinear temporal evolution equation of distribution is 
obtained that includes a self-consistent growth term. The derived equation is generally 
solved analytically by means of eigenfunction expansions. By focusing on the case with 
linear relaxation, two examples are considered as applications of the proposed general 
formalism. First, by assuming that the growth rate of a cell increases linearly with the 
state value x, the shift of the average state value x due to the growth effect is shown to 
be proportional to the variance of the state x and the relaxation time, similarly with 
the biological fluctuation- response relationship. Second, when there is a gap in the 
growth rate at some threshold value for the state x, existence of a critical gap value is 
demonstrated, beyond which the average growth rate starts to increase. This critical 
value is again obtained in terms of the relaxation time and the variance of x, all of 
which are experimentally measurable quantities. The relevance of the results to the 
analysis of biological data on the distribution of cell states, as obtained for example 
by flow cytometry, is discussed. 

Keywords: fluctuation, cell growth, distribution function, Fokker-Planck equation, flow- 
cytometry 
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1. Introduction 

Biological systems suffer fluctuations. No intracellular biochemical process can avoid 
fluctuations, because they arise from the motion and reaction of molecules. For example, 
gene expressions or abundance of some proteins in a cell fluctuate in time or by cells, even 
if they are measured at the same time after a cell division, for cells with identical genes 
(clones). Indeed, Elowitz has explicitly measured the numeric fluctuations of proteins in 
Escheria coli, by distinguishing intrinsic and extrinsic fluctuations ij. Such intracellular 
fluctuations have attracted both theoretical and experimental attention 01311100111, 
while the significance of the phenotypic fluctuations for adaptation [Hj and evolution 
[5J HO] has also been investigated. 

In general, let us consider the fluctuation of some quantity x characterizing the 
state of a cell, such as the number of proteins or gene expressions. Now, as a result of 
intra-cellular dynamics, x fluctuates among cells or in time. Let us denote the single-cell 
distribution of x by P S i ng i e (x). In principle, it can be obtained by repeating a single-cell 
measurement over an ensemble of cells. 

Here, however, we must be careful about the choice of the initial ensemble itself for 
such distributions. The initial distribution of cells chosen for an experiment depends 
on whether the cell can proliferate or not and the speed of cell replication, which may 
depend on the cell state x. Consider, for example, taking an ensemble of cells from a 
culture. Then the probability of choosing cells that have higher replication speeds will 
be larger, and the initial distribution of x will be biased accordingly. 

This problem is prominent in the measurement of cells from continuous cultures 
using flow cytometry or some other means [TT]. In flow cytometry, the characteristics 
of each cell (e.g., the magnitude of fluorescence when a fluorescent protein gene is 
introduced) are measured over a huge number of cells. It is now established as a 
standard, powerful tool to measure the distribution of states of cells. Here, if the growth 
rate of a cell is independent of the quantity x, the choice of cell ensemble is not biased 
by the value x, and thus the observed distribution P(x) by flow cytometry is simply that 
given by the distribution P S i ng i e {x). On the other hand, if the growth rate depends on 
x, the distribution P(x) may be altered from the distribution from single-cell dynamics. 

As an illustration, consider the case in which P S ingie{x) is a Gaussian distribution 
around x = Xo, while the replication rate of a cell increases strongly with x for x > Xq, 
assuming that x represents the abundance of some chemical that mediates the growth 
of the cell. In this case, it is naturally expected that the observed distribution P(x) 
should be biased towards x > xq. 

In general, the distribution P S i ng ie(x,t) has been studied with the use of stochastic 
processes to characterize the intra-cellular dynamics of the state x. Established 
mathematical tools such as Master's equation, Langeving's equation, and the Fokker- 
Planck equation [T2J EH are applied for such studies. On the other hand, as a biological 
unit (cell) replicates, the number of cells increases accordingly. This effect of replication, 
then, must be incorporated with these stochastic processes, to include both the single- 
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cell fluctuations and the growth dynamics of the cells together. 

Recently, there has been growing interest in exploring the relationship between the 
fluctuations of intracellular state and the response of the state to the change in external 
conditions, both theoretically and experimentally El El El- For example, a 
change in the concentration of some protein (or gene expression) against the change 
in the external condition (e.g., concentration of some chemical in the medium) may 
be measured experimentally, from which the response of such intracellular state to the 
environmental change must be unveiled. Here, however, the growth speed of a cell 
generally depends on the intracellular state, e.g., the abundance of such protein, because 
the protein is important to the function of the cell. Hence, the measured change of the 
protein concentration in response to external change involves both the internal change 
of the intracellular state and the change in the cell number distribution caused by the 
state-dependent growth rate. Thus, we should develop a theoretical tool to distinguish 
the two effects, based on the measurable quantities. In the present paper, by setting up 
an equation for P(x, t) that takes into account both the intra-cellular stochastic process 
and the state-dependent cell reproduction rate, we address this issue. 

We first derive the evolution equation of the distribution P(x, t) by extending the 
Fokker-Planck equation to incorporate state-dependent growth. (In the present paper, 
'growth' means the replication of a cell, and the replication rate in time is called the 
growth rate). The derived equation includes a term for state-dependent growth, from 
which is subtracted the average growth rate over all cells, leading to a source/sink 
term that corresponds to the growth/death process of a cell. The average growth rate 
gives a self-consistent term that is nonlinear in distribution P(x, t), but we can formally 
solve the equation through the eigenvalue properties of a Sturm-Liouville-type operator. 
After giving a general formulation of the equation, we present two simple examples of 
this formulation, by assuming the linear Langevin equation for the single-cell dynamics 
of the state variable. First, by considering the linear dependence of the growth speed 
on x, we obtain a formula for the shift of the average value of the state x. The shift is 
proportional to the product of the variance of the state, the relaxation time, and the 
proportion coefficient of the growth speed with x. For our second example, we study the 
case in which there is a threshold value of the state x for growth, and derive a formula 
for the change of P(x, t) to 'feel' the state-dependent growth. Cautious remarks are 
made on the interpretation of the distribution obtained from flow cytometry, while the 
relevance of our theory to evolution is also briefly discussed. 

Note that we do not discuss specific mechanism for the cell growth here. Rather, 
we introduce a function characterizing state-dependence growth generally and derive 
the distribution function. 
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2. Derivation of the equation for the distribution of cell state with 
reproduction 

Let us first introduce a variable x, which represents a state value of a cell, for example, a 
concentration of some chemical (or its deviation from the mean value). We assume that 
the temporal evolution of variable x in a single cell obeys some Markovian dynamics, 
that is, the value of x at time t is determined only by the value of x at some previous 
time. (Although biological systems may often retain some memory, this assumption can 
be acceptable as a first-step approximation, and indeed is adopted for most models.) 
Based on this assumption, we consider the following Langevin equation, which is often 
adopted: 

^ = -/(*(*)) + y/g(x(t)) V (t), (1) 

where / and g are functions of x(t), that govern the dynamic behavior of the variable 
x (the function g must be non- negative for all x); roughly speaking, the function / 
represents the force acting on the value toward its mean value and g represents the 
strength of the diffusion at the value. rj(t) is a Gaussian white noise term having the 
statistical properties: (r](t)) = for any t and (rj{ti)r](t 2 )) = 25 {t\ — t 2 ) for any t\ and 
t 2 . The distribution function P sing i e (x,t) indeed obeys the Fokker-Planck equation $ 
derived from the Langevin equation ((TJ) [T21 HBl ITTj 



d 



Psingle(% i t) ■ (2) 



dt dx 

We now introduce the growth (replication) of the cell, whose rate \i is dependent 
on the state value of x in the cell, and is a function of x(t), denoted by fj,(x(t)). To 
derive the equation for the distribution P(x, t) for this growth rate of the cell, we first 
write down the change in the distribution function at time t + At, given P(x, t) at time 
t, as: 

P(x,t + At)= W(x,x', At)P(x',t)(l + n{x')At)dx' , (3) 

J Xl 

where the term (1 + ^(x')At) indicates the effect of cell growth, while W(x, x' , At) 
is the transition probability that the system changes from the state with x' to that 
with x during the time interval At, which is determined by the Langevin equation (|T|) 
uniquely. Because of cell growth, the distribution function P(x, t + At) obtained above 
is not normalized in general, while the distribution P(x, t) must be normalized. To 
obtain the correct form of P(x,t + At), then, we must renormalize this distribution as: 

^dx'W{x,x',At)P{x',t){l + ^{x')At) 
C dxf% dx'W{x, x', At)P(x>, t)(l + ii{x')At) 

~(fx(x)-fl(t))P(x,t)At+ r W(x,x',At)P(x',t)dx', (4) 

J X\ 

I Here we have adopted Ito calculus; for Stratonovich calculus, one can simply replace / by / — 
in equation @. If g is constant, there is no difference. 
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where we have used the property of the transition probability, f^ W(x,x' , At)dx = 1 
for any x' and any At, and retained in the second line only the terms up to the first 
order in At. Here fx is defined by 



fi(x)P(x, t)dx, 



(5) 



which gives the mean growth rate of the cells at time t. In equation P|). taking the limit 
At — > and recalling the second term is reduced to the form of equation (J2J), we obtain 



dP(x,t) 



(n(x) - p{t))P(x,t) + 



d 



9 / \ 

+ d~x g{x) 



P(x,t). 



(6) 



dt dx 
This is the equation we desired to derive, the time evolution equation for the distribution 
function with x-dependent cell growth rate. 

As in the standard Fokker-Planck equation for the probability, we take the no-flux 
boundary condition as: 



d 

f( x ) + fa9{x) 



P(x,t) 



0. 



(?) 



If [i(x) = constant, i.e., for x-independent cell growth, the first term in equation 
(0), (p(x(t)) — p l (t))P(x,t), vanishes and accordingly equation © is reduced to just 
the usual Fokker-Planck equation (j2J); the influence of the state-dependent cell-growth 
appears only in the term (fj,{x(t)) — p,(t))P(x,t), which plays the role of source (sink) 
in the distribution density, if the growth rate at some point x is greater (smaller) than 
the mean growth rate, ft. Equation © obtained above is nonlinear in P because the 
term fi(t) involves P itself, so that it first looks rather difficult to analyze. Fortunately, 
however, the analysis turns out not to be so difficult, as will be shown in the next 
section. 



3. Analysis of the evolution equation of the distribution with growth 



In this section, we examine the structure of equation (jj 
and eigenvalues. We first introduce a linear operator 



» {X) + Yx 



d 

f( x ) + Qff9{x) 



with the aid of linear operators 



(8) 



and rewrite equation (0) as 
dP(x,t) 



dt 



-p(t)P(x,t) + L(x)P(x,t). 



(9) 



As the operator L is of the Sturm-Liouville type, we can, in principle, find all of 
its eigenvalues and corresponding eigenfunctions, and all the eigenvalues are real jTS] . 
We denote the eigenvalues and the corresponding eigenfunctions by Aj and 4>i{x), where 
the index i runs over the non- negative integers, i = 0,1,2,.., and the eigenvalues are 
ordered so that A, > Xj for i < j. From the definition, Aj and <pi{x) satisfy the relation 

L{x)<j>iix) = Xifatx). (10) 
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In general, we can introduce the adjoint operator of L, denoted by L\ and introduce 
the "left" eigenfunctions of L* denoted by ipi(x), for the eigenvalue Aj. As is well 
known, left and right eigenfunctions for different eigenvalues are orthogonal and can be 
normalized as J* 2 ipi(x)(f)j(x)dx = <5y, where is the Kronecker delta (Sy = for i ^ j 
and 5ij — 1 for % — j). 

With these relationships, we expand P(x,t) in terms of these right eigenfunctions 

as 

oo 

P(x,t) = Y / a 3 (t)<P 3 (x), (11) 

3=0 

where the {a,i(t)} are expansion coefficients that are related to the integration 

a 3 (t)= ipj(x)P(x,t)dx. (12) 

J Xl 

Next, we will express Jx{t) in terms of {a{\ and {Aj}. From the definition (JHJ) of ft, 

Ji(t) = / ji{x)P{x,t)dx — y^Ajaj(t) / <pi{x)dx 1 (13) 
Jxi i=0 Jxi 

where we have used the relation (JBJ), the boundary conditions (JZJ), and the relations 
(fTUj) and (JTTJ, successively. 

The time evolution equation for ai(t) is straightforwardly obtained by inserting (|TT| 
into (jnj), multiplying by ipi(x) and integrating it over x: 
dn ■ (t) 00 r x 2 

^1 = (A, - £A,^(t) / ^ds)^*) (14) 
at j_ Jx\ 

In summary, the partial differential equation (JHJ) for P is reduced to a set of ordinary 
differential equations for {aj}, while the initial conditions of a, are given from the relation 
(JHJ): aj(t ) = J^ 2 vpi(x)P(x,t )dx for the initial time t - 

Note that there remains a freedom in the choice of <pi(x) and if>i(x), because 
the normalization condition is still satisfied under the change of 4>i(x) — > Ci4>i(x) and 
ipi(x) — > (l/ci)if>i(x) with any constant q ^ 0. By taking advantage of this freedom, we 
can introduce, for convenience, another normalization condition: 

/ <pi{x)dx = 1 (15) 

J Xl 

for all the right eigenfunctions whose integral over x does not vanish. Indeed, 
this normalization is easily achieved by re-scaling the eigenfunctions ipi(x) — > 
ipi(x) f^ <j>i(x')dx' and 4>i( x ) ~* (pi{ x )/ fxl (i>i{ x ')dx' . If /j? 2 4>i(x')dx' vanishes, we simply 
leave the original eigenfunctions, and we call eigenfunctions with J^ 2 (pi(x')dx f = 
"non-contributing eigenfunctions". Note that for the 0th right eigenfunction, O > this 
normalization is always possible, because the 0th right eigenfunction does not take 
4>o(x) = for any x JB]- With this choice of normalization, equation (|T4*|) is simplified 
as 

^ = (A,-^V 3 (M), (16) 



On the Distribution of State Values of Reproducing Cells 



7 



where the prime over the summation symbol indicates that the summation is taken over 
all eigenfunctions except non-contributing ones. 

Equation (fTHjl tells us that any eigenfunction (j>i{x) of the linear operator L, except 
for the non-contributing ones, gives a stationary solution of equation (JHJ), because any 
set { a,i(t) = 1 and aj(t) = for j ^ i } is a stationary solution of (JHJ). Among those 
stationary solutions, however, only the solution with a,j(t) = 5j t0 , is stable. 

To show this, we make a linear stability analysis of these solutions. Consider the 
solution a,i(t) = 5ik for given k, and introduce a perturbation 5ai(t) as a,i(t) = 5ik + Sai(t) 
(i = 0, 1, ...,). Then, inserting this into and retaining only the terms of first order 
in Sa, we obtain 



The eigenvalues of the matrix {A^} are easily shown to be (Aq — A&), ... , — A&, 
(Afc+i — Afc), ... Recalling that the eigenvalues are ordered so that A* > \j for i < j, we 
can easily show that all the stationary solutions for k > are unstable, while if Ao > 
the solution with k = (i.e., with etj = 5 i0 ) is stable. In other words, only the mode 
with the largest growth rate remains as a stationary solution, as is expected. 

The requirement Ao > for the stability of the system is quite reasonable. 
Otherwise, all Aj are negative, which means there is no growth at any state, and all 
the cells would become extinct with time (recall that Aj is equal to the growth rate 
of the mode represented by the ith eigenfunction). To have a positive growth rate for 
the stationary distribution, Ao > is therefore necessary. The condition Ao > simply 
means that the cells (or units) continue reproduction without extinction. 

Now, the stationary solution of equation (JHJ) is given by 4>o{x), the eigenfunction 
of the operator L corresponding to the maximal eigenvalue Ao- Similarly to the case of 
the standard Fokker-Planck equation, the eigenvalue problem of the operator L can be 
transformed into that for the Schrodinger-type equation whose "potential" is given by 
the functions f(x), g(x), and fj,(x) (see Appendix A). Hence we can use the methods 
and solutions developed in quantum mechanics. 

4. Two simple examples of the evolution of the distribution 

In this section we study two simple examples of equation (JHJ) by linear or threshold-type 
dependence of the growth rate on x. We choose f(x) = kx and g(x) = D in equation (JHJ) 
with k and D positive constants; the reasons for this choice are: (i) that the Gaussian 
distribution is often observed to be the stationary distribution of a biological state, 
while this linear Langevin equation is the simplest to realize the Gaussian distribution 
(the log-normal distribution is sometimes observed in cells El HHj , but in this case we 
can simply use the logarithm of the quantity as the variable x that concerns us), and 
(ii) that this linear Langevin equation has been thoroughly investigated in physics and 



dSai(t) 



oo ' oo ' 



dt 



(Aj - X k )Sai(t) - 5 ik Y^ ^j3aj(t) = Aij8aj(t). 

j=0 j=0 
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mathematics; it models the motion of a Brownian particle in a harmonic potential, so 
that we can easily see the effect of the state-dependent growth introduced here. 

4-1- fi{x) linearly dependent on x 

We study the case (i(x) = ax + b for x to [—00, 00] in equation ©, where a and b are 
constants. It is natural to study the linear case as the simplest non-trivial example. 
Indeed, as long as the range of x in concern is small, gradual change in u(x) can be 
approximated by linear change. 

In this case, we can obtain all eigenvalues and their corresponding eigenfunctions 
of L as A n = + b - kn and <f> n (x) = N n H n (^(x - exp[-^(x - - 
j^x 2 ], where H n (x) is the nth Hermite polynomial in x and N n is the normalization 
constant determined by the normalization condition ()15p. In particular, the stationary 
distribution is obtained directly as 

o (x)=iV o exp[-A (x _^) 2] , (17) 

while the temporal evolution of the distribution is obtained with these eigenvalues and 
eigenfunctions and with the reduced equations (fTHjl for {a^}. 

Fortunately, however, in this case there is a more convenient way to obtain the 
dynamics of the system: if the system starts with a Gaussian distribution at some 

initial time, the temporal evolution of the distribution preserves the Gaussian form. By 

(g-g(Q) 2 

taking a Gaussian distribution P(x, t) = , e 2 /*w with a and 8 as the mean 

value and the variance, it can be shown (see Appendix B), that the temporal evolution 
preserves the Gaussian form when the time evolution equations for a and B are given 
by ^ = adit) - ka(t) and * = -2k0{t) + 2D. 

These equations indicate that while the temporal evolution of the variance is 
completely the same as the case for a constant u, the evolution of the mean value 
is influenced by the state-dependent growth; the mean value is shifted in the direction 
of larger u, driven by its variance. In the stationary state, as is also given in equation 
(|17|). the mean value (peak position) shifts with the degree aD/k 2 compared with the 
case without the growth term (or, from the case with constant « (i.e., a = 0)). Note 
that this change in the mean value in the stationary state is proportional to the variance 
of the original distribution, which is given by D/k, i.e., 

aD a , 9 \ 



k 2 k 

where (...) is the average of the stationary distribution P(x), and 5x = x — (x). 

In other words, the larger the variance of the distribution is, the more the mean 
value shifts. Correspondence with the fluctuation-response relationship [201 |H] is 
interesting, because the shift in the growth is proportional to the original fluctuation. 
In addition, response to a higher growth state is possible only under the fluctuation of 
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the state, which demonstrates the relevance of phenotypic fluctuation to adaptation. 
With this shift of Ax, the average growth rate of a cell changes with 



which is an experimentally measurable quantity. Hence, the right hand side of equation 
(JTHJ) is represented by measurable quantities, because k is simply the relaxation time, a 
is estimated from equation (|19|) and the variance {(Sx) 2 ) is measurable. 

4-2. A threshold for growth: the step function fj,(x) 

We consider equation © with fj,(x) = a, Q(x — xo) +b, where a, b, and xo are constants, 
and is the so-called Heaviside step function; 0(x) = for x < and Q(x) = 1 for 
x > 0. We study this case, because in biological systems, a threshold for reproduction 
sometimes exists. 

In this case, the eigenfunctions are written analytically with the use of confluent 
geometric series and the corresponding eigenvalues are obtained, by transforming the 
equation to the Schrodinger equation (see Appendix A). Because the complete analytic 
form is rather complicated, we discuss only the results of numerical calculations here. 

First, we consider the stationary distribution of equation When the position xq 
of the step of fi(x) is within the standard deviation of P S i ng ie(x), i.e., < x < (we 
consider only the case of non-negative xq), the stationary distribution gradually moves 
toward the position Xq, as the parameter a increases. On the other hand, when the 
position xq is outside the standard deviation of P S ingie(%), i-e., x > J^, the stationary 
distribution does not change much until the parameter a reaches some critical value 
a c . As a increases beyond that value, the distribution shifts smoothly to larger x. The 
existence of the critical value a c is demonstrated in figure which is a plot of the 
total amount of the distribution in the right region (x > Xq) against the relative growth 
rate a (see figure (J2J)). 



The critical value of a c is estimated to be a c ~ kxoy |j, as is confirmed numerically 
(see inset of Fig Indeed, this value of a c coincides with the inverse of some 

characteristic time, that is the average time required for a cell in a higher-growth 
state (x > x ) to change to the lower-growth state (x < x ). This numerical result 
is reasonable: if the relative growth rate a is smaller than a c , cells change to the state 
x < Xq before they grow sufficiently in the higher-growth region x > xq. The cells 
cannot 'feel' the higher-growth region, so that the difference in growth rates does not 
influence the cell population distribution. 

Next, we briefly explain the dynamic behavior of the distribution when the relative 
growth rate is greater than a c and the distribution is initially localized at x < xq. To be 
specific, we set P(x,t ) = 5(x), i.e., localized at x — 0. The temporal evolution of the 
distribution is given in figure Q. Here: (i) first, the distribution behaves as if it does 
not 'feel' the state-dependence of fi(x), until its tail touches xo, the edge of the step 
function, (ii) After the tail of the distribution reaches the edge of the step function, the 
distribution in this tail region starts to grow faster (see figure (JH1)); at this stage, the 



A/i = aAx, 



(19) 
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distribution has two peaks, (iii) Finally, the distribution converges to a single peak at 
the mean value at around Xq, the position of the step of fi(x). This temporal evolution 
to a higher growth state is in contrast to the linear case, where a single-peak distribution 
is preserved and only the peak position is shifted. 

In the present example, the stationary distribution has a single peak. For some 
forms of f(x), however, the stationary distribution has two peaks, even though the 
single cell distribution (without the x dependence of [i(x)) has a single peak. For 
example, for f(x) = 2sgn(x) with the present form of two peaks coexist (see figure 
(jU). Here, for large x (x > x (— 4)), the growth rate is high and the distribution is 
confined within some range, so that the distribution has one peak in that region, while 
for small x (x < Xq), not all cells grow so that the distribution of the cells tends to 
decrease. However, many cells that have grown in the higher- growth region flow into 
the lower-growth region because of the effect of the force of /, so that the distribution 
has another peak there. 

5. Conclusion and discussion 

In the present paper we have posed the question of how the distribution of an 
intracellular state variable (say the abundances of some chemical or degree of gene 
expression) is altered due to the state dependence of the replication rate of a cell. 
To discuss the temporal evolution of the distribution of the internal state x of 
such replication units, we have incorporated the state-dependent growth rate into 
the standard Fokker-Planck equation. By considering the population distribution of 
replication units with Langevin equation dynamics, we have derived a general equation 
for the temporal evolution of the distribution of states P(x,t). The derived equation 
includes a self-consistent term arising from the growth rate. In spite of this non-linear 
term, we can formally solve the equation as an eigenvalue problem of the Sturm-Liouville 
type. Note that the formalism presented here is rather general, as is the Fokker-Planck 
equation. 

After giving a general analysis of the equation, we have studied two simple examples, 
assuming the linear Langevin equation for single-cellular dynamics. First, when the 
growth rate increases linearly with the state value x, the average of x over cells increases 
in proportion to its variance, which reminds us of the fluctuation-response relationship 
in physics, while the proportion coefficient is estimated by the increase of the growth 
rate and the relaxation time. Note that the shift of population distribution to a higher 
growth state is possible only with the fluctuation of the internal state. Our result implies 
that the response of x to environmental change is proportional to its variance. In other 
words, fluctuations in chemical concentration, which have been studied extensively, are 
relevant to biological adaptation. 

Now let us return to the question raised in the introduction. We measure an 
intracellular state variable (x), from an ensemble of cells, and study its change against 
the change in external conditions. Here we change the environmental condition (e.g., 
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nutrient concentration) and the cell state value x (e.g., the concentration of some 
enzyme) is changed accordingly. After the cell distribution becomes stationary, we can 
measure this change of the average x denoted by (Ax) total that is caused by the change 
in the environmental condition. Now, from this measurement, we are often interested 
in detecting the change in the stationary state of x, to explore intracellular dynamics. 
However, such an intracellular state variable x is often also related with the ability 
for cell growth. Hence (Ax) total is also influenced by the change in cell growth speed, 
and this may deviate from the change caused by the intracellular dynamics (Ax) single . 
Then, can we estimate the change of the internal state (Ax) single from the measurement 
of (Ax) total 7 If we confine our discussion only to the linear regime, we find 



from equation (fTHj) . Here the latter term can be estimated from the standard 
measurements. First, through equation ()19|). a can be estimated from the change in 
the average growth rate of cells. Second, k is simply the relaxation time. Hence, by 
measuring the temporal change of (x(t)), and by fitting the approach to its stationary 
value by an exponential form, one can estimate k. Finally, from the variance of the 
state value x at a stationary state (by flow cytometry or by other means), we can obtain 
((5x) 2 ). Accordingly, we can estimate the term | {(Sx) 2 ), so that the intracellular change 
of x is estimated from the observable quantity (Ax) tot . 

In our second example, we studied the case with a threshold-type dependence of 
the growth rate on the state x. When the position x of the step of ^(x) is outside the 
standard deviation of P S in g ie(x), i.e., when x > y^f, the distribution does not change 
significantly until the relative growth rate a reaches a critical value a c , beyond which the 
distribution starts to shift to the higher-growth region. From the biophysical viewpoint, 
the value a c corresponds to the inverse of the average time required for a cell to change 
from the higher-growth state (x > Xq) to the lower-growth state(x < Xq). 

Here we have found that the distribution of the state variable often exhibits double 
peaks over a long transient time. For some form of f(x) and fi(x), a double-peak 
stationary distribution is also obtained, even if P S i ng ie(x) has only a single peak. This 
raises a cautious remark on the interpretation of the distribution observed in flow 
cytometry. Even if double peaks are observed, this does not necessarily mean that the 
internal cell dynamics (e.g., gene expression network dynamics or metabolic dynamics) 
have bistable states. One of the peaks may be associated with the flow of population 
due to the difference in reproduction speeds. 

Several extensions of the present formulation are straightforward. Although we 
mainly discussed the case with a single state variable, extension to a higher-dimensional 
case is straightforward. Inclusion of a memory term to go beyond Markovian dynamics 
will be possible, although we expect that most of the results on the linear and step- 
function cases above are still valid in the non-Markovian case. 

Although we have given our formulation here for a reproducing cell with an internal 
state (e.g., chemical concentration), the present formulation can be applied generally 




(20) 
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to any reproducing system with a growth rate dependent on its internal state. For 
example, it can be applied to an artificial cell or a replicating biochemical system with 
a growth rate that depends on its internal catalytic activity. Furthermore, application 
to continuous evolution is possible. By taking x as a Hamming distance from a typical 
gene, the evolution process to change x to a given phenotype having some function 
can be considered. Here, the reproduction rate depends on x, which gives fi(x), while 
the diffusion process in x is simply the mutation, with D as the mutation rate. As 
non-functional mutants are more common, the mutation in the change of function (or 
activity) has a drift to a smaller regime, leading to a 'force' term towards x = as in 
equation^). The temporal evolution of the distribution of gene x is thus analyzed by 
using our equation (JHJ), while in some examples, the steady state with positive growth 
rate collapses ^H], with the increase of the mutation rate, as the largest growth speed 
Ao becomes negative, which leads to error catastrophe. 

A biological unit reproduces at a rate that depends on its state. The present Fokker- 
Planck equation with growth and death provides a basic equation for such problems in 
general. 
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Appendix A. Transformation of the linear operator L to a Hermite operator 



In this section we transform equation (JHJ) to a type of Schrodinger equation, to show 
explicitly that the operator L defined by (JBJ) is transformed to an Hermite operator. 
Here we follow the standard transformation from the Fokker-Planck equation to the 
Schrodinger equation except for the existence of the terms concerning fi(x). 

We first introduce a new variable y defined as y(x) = f£ Q yJ^Zdx', where x is some 
number on [rci,^]. According to this transformation, the distribution can change to 
P(y,t) = QrfteP(x,t) = yfa$P(x,t). With these new variables, we can write equation 
© as: 



P(y,t) = -jl(t)P(y,t) + 



Ky) + 



d_ 

dy 



f(y) + D 



d_ 

dy 



P(y,t), 



(a.i; 



where f(y) = y (/fa) + 2#'fa)) anc ^ A(l/) = M fa (?/))• 9'( x ) 1S the derivative of g with 
respect to x and x(y) is the inverse of the function y(x). Note that Jjl does not change 
by this transformation. 

By further introducing two new quantities $(?/) 
P(y,t), equation (|A.1|) is rewritten as 



jyltfdy' and tf(y,t) 



»(«) 
e 2 



dt 



p,(t)V(y,t) + 



V(y) + D 



dy 1 



fi(t)y(y,t)+H(y)V(y,t) 



(A.2) 
(A.3) 
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where V(y) = fi(y) 



£^ + m andH{y) = [ v{y) + D *_ 



The operator H obtained 



above is evidently a Hermite operator, and indeed the eigenvalue problem of H^f is 
simply a type of Schrodinger equation. Accordingly, the exact solutions or techniques 
developed for Schrodinger equations can be applied to our problem. 



Appendix B. Temporal evolution preserving a Gaussian distribution for the 
linear fi(x) case 

When f(x) = kx, g{x) = D, and = ax + b, equation (JHJ) becomes 

d 



dP(x, t) , i \ \ 7 — / \ d 
- a(x - (x) t )P(x,t) + 



kx + D 

ox 



P(x,t), (B.l) 



dt dx 

where we have used the normalization condition P(x,t)dx = 1, and have adopted 
the notation (...) t = f^ ...P(x,t)dx. Multiplying both sides of equation (|B.1|) by x and 
x 2 and integrating each case over x, we obtain 

^Mi^t-ix} 2 )-^ (B.2) 

^M(x 3 ) t - (x) t (x 2 ) t ) - 2k(x 2 ) t + 2D. (B.3) 
Suppose now that the solution of equation (|B.1J) is a Gaussian distribution, i.e., 



1 (x-q(t)) 2 

P(x,t) = —. e urn , (B.4) 

where a and (3 correspond to the mean value of x and its variance, respectively, which 
are related to (x) t and {x 2 ) t as a(t) = (x) t and (3{t) = (x 2 ) t — (x) 2 . Using equations (|B.2|) 
and ()B.3|) and the property of the Gaussian distribution (x 3 ) t = 3a(t)j3(t) 2 + 3a(t) 3 , we 
can derive the time evolution equation of a and f3 as follows: 

^La/3(t) - ka(t) (B.5) 
dt 

d(3(t) d(x 2 ) t d(x) t 

dt dt { H dt 

=- 2k(3(t) + 2D. (B.6) 



On the other hand, inserting the form of (jB.4|) into equation (jB.l|) and simplifying 
the equation, we obtain the equation 

2{x -a(t))f3(t)(ka(t) - a(3{t) + + ((x - a(t)) 2 - P(t))(-2D + 2k(3(t) + = 0. 

(Jib (U/L 



The time evolution equations of a and j3 satisfy the above equations (jB.5|) and (jB.6J) . 
and the Gaussian distribution is the solution of equation (jB.l|) (as the solution with 
temporal evolution is unique). 
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Figure 1. The total amount «2 of the distribution in the higher-growth region 
(x > xq — 3), calculated as n-i — P(x)dx, against the relative growth rate 
a. This plot clearly indicates the existence of the critical value of a, a c , which is 
defined here as the value where ri2 = 1/10 . The inset, a plot of the critical value 
of a c against the position Xq of the step of n(x), shows the dependence of a c on 

xq, which is fitted well by: a c ~ kxoy^. These calculations were carried out for 
f(x) — x, g{x) — 1, fi(x) — aO(x — 3), for the range of x, [—6, 6]. 




Figure 2. Some profiles of the stationary distributions for different relative growth 
rates a — 2 (black), 5 (red), and 8 (blue). We can see that the distribution for a = 2 
is hardly influenced by (J,(x). We choose the same equation as for Figure 1, i.e., 
f(x) = x, g{x) = 1, fx(x) = aQ(x — 3), with the range of x, [—6, 6]. 




Figure 3. Temporal evolution of the distribution for equation JBJ for /i(x) = 
20 0(2; — 3), g(x) = 1, f(x) = x, with the range of x, [—6, 6]. The initial condition 
is given by P(x,to = 0) = 5(x). The black, red, and blue curves show the 
distributions at t = 0.6, 1.04, and 1.6, respectively. The double-peak distribution 
is observed during an intermediate period. 
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Figure 4. An example of the stationary distribution of equation © having double 
peaks, for f(x) = 2 sgn(a;), g(x) = 1, (i(x) = 2.4 0(x — 4), with the range of x, 
h7,7]. 



